{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Data Analysis, 3rd ed\n",
    "## Chapter 2, demo 1\n",
    "\n",
    "Authors:\n",
    "- Aki Vehtari <aki.vehtari@aalto.fi>\n",
    "- Tuomas Sivula <tuomas.sivula@aalto.fi>\n",
    "\n",
    "Probability of a girl birth given placenta previa (BDA3 p. 37).\n",
    "437 girls and 543 boys have been observed. Calculate and plot the posterior distribution of the proportion of girls $\\theta $, using\n",
    "uniform prior on $\\theta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# import necessary packages\n",
    "\n",
    "# plotting\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# fast array routines for numerical calculation\n",
    "import numpy as np\n",
    "# scipy contains various scietific tools, such as beta distribution\n",
    "from scipy.stats import beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# add utilities directory to path\n",
    "import os, sys\n",
    "util_path = os.path.abspath(os.path.join(os.path.pardir, 'utilities_and_data'))\n",
    "if util_path not in sys.path and os.path.exists(util_path):\n",
    "    sys.path.insert(0, util_path)\n",
    "\n",
    "# import from utilities\n",
    "import plot_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# edit default plot settings\n",
    "plt.rc('font', size=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The posterior distribution is Beta(438, 544). Plot the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Uniform prior -> Posterior is Beta(438,544)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEGCAYAAACHGfl5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VOXd//H3mZnsk30hCSSBsO9LAFFAwSpWUHyw/QmK\n8GjVll61rq1WtAIWa921RR93qAqurSsqomwi+yYBAiEkkJAQspF9ncz9+2OSMQkJ2XNm+b6uK5o5\n63dODvOZc+5z7qMppRBCCOHeDHoXIIQQQn8SBkIIISQMhBBCSBgIIYRAwkAIIQQSBkIIIZAwEEII\ngYSBEEIIJAyEEEIApnZOL7crC+GIVs6y/f/WtfrWIVqi6V1Aa+TIQAghhISBEEIICQMhhBBIGAgh\nhEDCQAghBBIGQgghkDAQQgiBhIEQQgjaf9OZED2q2mIl+WwJiZlFHD9byozhvZgUH6p3WUK4HAkD\n4ZBO5JbywMcHSTxdSHWt7cZ3owartqVx1+UD+OMvBmE0OPxNnUI4DQkD4XByS6q45a1dlFbWMHdc\nJEN7+TG0lx8hfh78Y30aL3yfwq60Av550zjCzF56lyuES9CUald3Q9I3kehW5dUWbnxtB8eyS3h1\n7lCGR5kbjVdK8VliLk9vOEmAtwcvzx/HhH5y2kj6JnJ4Dn8YKw3IwmHUWhV3vbefxMwiHr9mwHlB\nAKBpGv8zKoJV80fgZdL43Tt7KSqv1qFaIVyLhIFwCEopln5+mO+ScvjT5XFcNiD4gtMPDPfliWsG\ncq68hqe/SeqhKoVwXRIGTuqhhx7ihRdeaNO0S5cuZenSpd1az8SJEzl8+HCH5/9oz2ne2XGKBROi\nuGFsZJvmGRrpx6/HRLBm92kOnT7X4XULISQMnFJubi5vv/02v/vd7xoN379/P5MnT8bX15eJEyeS\nnp7eYzX96U9/4tFHH+3QvFWWWp7/LpnR0Wb+eGlMu+b9/ZQYArxNPPLpIaxWadISoqMkDJzQqlWr\nmDlzJj4+PvZhp0+fZubMmTz44IPk5+cTHx/P8uXLe6ym2bNns3HjRrKzs9s970d7TnOmqJI7LumD\nQWtfO1uAt4m7Lo3lwOliPtx9qt3rFkLYSBg4qNWrV3PxxRczd+5coqKiiImJ4euvvwbg66+/5rLL\nLms0/f33388dd9zB7Nmz8fHxYd68eezevfu85ZaWlmI0Gjlz5ox92KFDh4iKiqKkpKTFel599VVm\nzpzJH/7wB8LCwoiOjmb9+vX28d7e3iQkJLBu3bp2vc9qi5WXN6UwKtrMRXEB7Zq33jUjwhgVbebJ\ndccolMZkITpEwsBBJSYmcuDAAebOnUtmZiZ33303ixYtso8bPHiwfdri4mI+++wzbr/9dvswq9WK\nt7f3ecs1m80MGTKEffv22Yf95S9/YfHixfj7+7dYz08//cSOHTuYPXs2OTk5/O53v+PJJ59sNM3Q\noUP56aef2vU+P957mqzCSu64uDdaO48K6hk0jQev6EtRhYV/fNXxdgsh3JmEgYNKTEzk3nvv5frr\nr8dgMLBw4ULS09OprKyksLCw0Qf3999/T01NDaNGjSIoKIigoCDmz59PXFxcs8ueMGGCPQy2bNnC\nkSNHzmt/aOrgwYP85S9/4aqrrsJgMDBs2LDzpvH396ewsLDN77HaYuWljSmMiDIzqW9gm+drzuAI\nP341uhcf7c0is6CsU8sSwh1JGDioxMREfv3rX9tf5+TkYDab8fb2Jjg4uNEpnZMnTzJ79mwKCwvt\nP9OnT+eXv/xls8tuGAYPPPAAf/vb3/D09GyxFqUUiYmJXHvttfZhhw4dOi8QSkpKCAoKavN7/GT/\naTILK7i9E0cFDc0fH0mtgjU7T3Z6WUK4GwkDB1RYWEhGRgbh4eH2YR9//DFXX301AKNGjSI5Odk+\nrqqqCl9fX/vrtLQ09uzZw+zZs5tdfn0Y/Oc//6GyspKbbrrpgvWcPHkSi8XS6NTU/v37GTNmTKPp\nkpKSGD16dJveY02tlRUbUhgW6cfkfp07KqjXJ8ibi/sG8vG+LGostV2yTCHchYSBA0pMTMRoNLJm\nzRosFgtr167l5Zdftt8rMHPmTDZv3myffsKECWzevJmsrCwyMjK46aabePzxxwkJCWl2+aNHjyY7\nO5v777+fJ554otG38ltuuYVbbrml0fQHDx5k5MiRGAw/7y779+9v9MFfWVnJ3r17ufLKK9v0Hj/Z\nn0nGuYpOtRU05/rREZwtqWb94TOtTyyEsJOO6hxQYmIi8+fPZ/v27QQHBzN48GA+/fRT+2mZhQsX\nMmbMGCoqKvDx8eHyyy/nmmuuYdCgQYSGhvLggw9yxx13tLh8Ly8vRo4cidlsth9t1MvIyGDevHmN\nhh08eLDRUUBeXh7Z2dmMGDHCPuyLL75g2rRpREdHt+k9vrU1jcERvkyJb/tppbaYGh9EuNmD1TtP\nMXN0ny5dthCuTMLAASUmJjJmzBjuvffeZseHhYWxcOFCXn31Ve655x40TeOVV17hlVdeadPyq6ur\nycnJYcWKFecNz8rKOu/I4K9//et566+pqWk07JlnnuHNN99s0/qTzhRzNLuEB34R16VHBQAmo4H/\nGRnBG9szOZlbQt/wlq+QEkL8TE4TOaDExESGDh16wWn+/ve/c88993Ro+cuWLWPy5MlMmjSp0XBP\nT0+SkpLw8PBo9zJ37tzZ6EjhQj49kInJoHHl4O7pbfS6keFoGry7Pa1bli+EK5IwcECHDh1iyJAh\nXba8adOmMW3aNPbt20dgYCBbtmzhX//6V5ctvz2sVsVn+7O4uG8gwb7tD522iAzwYkp8EP89cIaq\nGmlIFqIt5DSRA2rPtfptMW3aNPvvRUVFXbrs9tqRmk92cSV3X9q95/N/NboXW04U8tVPp5kzvvn7\nLYQQP5MjA9GjPtmfiZ+nkan9L9xFdWdN6htIVIAna3b1XGd9QjgzCQPRYypravn60BkuHxSCt0f3\n7npGg8acURHsTi8mJVvfoyEhnIGEgegx3yWdpbSqlpnDeuYxlTOHhQHw5U+ZPbI+IZyZhIHoMZ/u\nzyTC7Mm4Ph3rnbS9IgO8GBbpx/qkHNr5rG8h3I6EgegR+aVVbDqWy1VDQzEaeu7Z4JcPDOFwdhkZ\n+aU9tk4hnJGEgegRaxPPYLEqZtWduukp0wfaGqq/Tszq0fUK4WwkDESP+GRfJgPDfRkQ7tv6xF0o\nLsSH+FAfvj2S06PrFcLZSBiIbpddVMn+jEJmDGm+47zudvnAEPafLianuEKX9QvhDCQMRLfbeMz2\nrfzSbr63oCXTBwVjVbAuUa4qEqIlEgai232flENUgBfxoT66rH9QuC+9A71Yd/isLusXwhlIGIhu\nVVlTy48peUyJD+ryHkrbStM0pg8MZuepIgrLKnWpQQhHJ2EgutX21HwqamqZ2r9rn1vQXtMHhlBT\nq+ShN0K0QMJAdKsNSTl4exhIiOmZG81aMjLaTKifh5wqEqIFEgai2yil2HA0h4mxAXiZ9N3VDJrG\ntAHBbD1xjvKqmtZnEMLNSBiIbpN8tpTMwopu76G0raYPDKHSYmXDETlVJERTEgai23x/1HZKZnI/\nfdsL6o2P8cfP02i/1FUI8TMJA9FtNiTlMKSXHxH+nnqXAtiejzwhNoBtqeewWq16lyOEQ5EwEN3i\nXFk1+9LPMcVBjgrqTeobyJniapLlGQdCNCJhILrF5uRcrAqm6HxJaVMX9w0EYNNRuapIiIYkDES3\n+P5oDiG+HgyL9NO7lEZ6B3kTE+TF1pR8vUsRwqFIGIguV1NrZfOxHCbHB2HQ6a7jC5nUN4g96cVU\nVMslpkLUkzAQXe5ARiHFlRamxDvWKaJ6F/cLpNJiZWdKrt6lCOEwJAxEl9t6PA+DBhNi9b3ruCXj\nYwIwGTQ2J0sYCFFPwkB0uR9P5DG0lx8B3ia9S2mWr6eRUdFmfjxRIM9GFqKOhIHoUmVVFg6kFzrs\nUUG9i/sFkpxbTnZhud6lCOEQJAxEl9qVVoDFqpgQF6h3KRd0cV9be8bmY3KJqRAgYSC62NaUPDyN\nGqOj/fUu5YIGRfgS7GNii7QbCAFIGIgu9mNKHqN7++Pt4di7lkHTmNQ3kO1phdTWStcUQjj2v1jh\nVPJKqziaXcLEOMduL6h3Ud9AzlVYOJhRoHcpQuhOwkB0mW0nbHf1Toh17PaCepPqu6aQdgMhJAxE\n19mWkofZy8jQXo7VBUVLwvw8GRTuy48n5MhACAkD0WW2puQxPiYAo8HxuqBoyfjYAA5mllBeWa13\nKULoSsJAdIn0/HJOn6tgooNfUtrU+NgAqmsVu9Ok4zrh3iQMRJfYmpIHwEQHv9msqXF9/DFosO1E\nnt6lCKErCQPRJX48kUe42ZO4EG+9S2kXs5eJIb382Jl2TrqmEG5NwkB0mtWq2J6Sx8TYADQH7LK6\nNeNjAjh0ppSySunSWrgvCQPRaUnZxRSU1zhde0G9CbEBWKyKnalyN7JwXxIGotO22+8vcK72gnqj\ne/tjNGhsS5F2A+G+JAxEp+1IzSc22JsIf0+9S+kQX08jwyP92HmyUNoNhNuSMBCdUmtV7EwrICHG\nsTuma82E2ACSskspLKvSuxQhdCFhIDol6UwxJZUWxsU4Z3tBvfGxAdQq2HFC2g2Ee5IwEJ2yI9XW\nXjCuj1nnSjpnZJQ/HkbN3v4hhLuRMBCdsiO1gJggb3r5e+ldSqd4exgYFW1m50m530C4JwkD0WFW\nq2L3yQLGOXl7Qb3xMQEk55STX1KhdylC9DgJA9FhSdnFFFXUMD7GOS8pbWpCbCAK2JYi7QaOZPXq\n1cyYMaPLl/vDDz8wePDgLl+uo9E07aSmaVe0Np2EgeiwHam2rp/H9nGNI4PhUX54mwzSbqCjkydP\nomkaFovFPmz+/Pl8++23Xb6uqVOncuzYsS5frrOSMBAdtjM1nz5BXkQGOHd7QT0Po4HRvc3sOiX3\nG9Rr+KHsSutyFZqmmbpqWRIGokOsVsWutALG9XGNU0T1xscGcCKvgtxi12036Nu3L0888QTDhg0j\nODiYW2+9lcrKSgA2bdpEnz59ePLJJ4mMjOTWW28F4PXXX2fAgAGEhIQwe/ZssrKy7MvTNI1//vOf\nxMfHExYWxp///GesVttzpa1WK8uXLycuLo6IiAgWLlxIUVER8PNRwJtvvklsbCyXX345l156KQBB\nQUGYzWa2b9/OqlWrmDJlin1927ZtY8KECQQGBjJhwgS2bdtmHzdt2jT++te/MnnyZPz9/ZkxYwZ5\nec3fWV7/Xhtul2eeeYZRo0YRGBjI3Llz7dulqdraWu6//37CwsLo168fK1asaHREU1RUxG233UZU\nVBS9e/dG07TlmqYZ67bXLZqmbdU07RlN085pmpamadrVDbZnoKZpb2qadkbTtMxm5v1R07TnNU3L\nB5ZqmtZf07QNmqbla5qWp2naak3TgtqwKzQiYSA65NjZEgorakhwkfaCevXhtt3F2w1Wr17NunXr\nOHHiBMnJySxfvtw+Ljs7m4KCAk6dOsVrr73Ghg0beOihh/jwww85c+YMcXFxzJs3r9HyPvnkE/bs\n2cO+ffv47LPPeOuttwBYtWoVq1atYuPGjaSmplJaWsqdd97ZaN7NmzeTlJTEunXr2LJlCwCFhYWU\nlpZy8cUXN5q2oKCAWbNmcdddd5Gfn899993HrFmzyM//+dTemjVrWLlyJTk5OVRXV/PMM8+0ebt8\n+OGHfPPNN6SlpXHw4EFWrVrV7HSvv/46X3/9NQcOHGDfvn18+umnjcbfcsstmEwmUlJS2L9/P8AM\n4PYGk1wEHAPCgKeAN7Wfe3lcBViAAcDYFuZNBXoBjwMa8AQQDQwFYoClbX7TdSQMRIf8fH+Ba7QX\n1BsWaWs3qH9/rurOO+8kJiaGkJAQHn74Yd577z37OIPBwLJly/Dy8sLHx4fVq1fzm9/8hnHjxuHl\n5cUTTzzB9u3bOXnypH2eBx98kJCQEGJjY7nnnnvsy1u9ejX33Xcf8fHxmM1mnnjiCd5///1Gp4SW\nLl2Kn58fPj4+rda9du1aBg4cyIIFCzCZTNx4440MGTKEL774wj7NrbfeyqBBg/Dx8eGGG27gwIED\nbd4ud911F9HR0YSEhHDttde2OO+HH37I3XffTZ8+fQgODuYvf/mLfdzZs2f56quveOGFF/Dz8yMi\nIgLgeaBhgp5SSr2ulKoF/g1EAb00TesFzATuUUqVKaVympk3Syn1L6WURSlVoZRKUUqtV0pVKaVy\ngeeAy9r8put02fkm4V52pOYTHehFVKBrtBfU8zAaGNXbzJ70IpRSTtkld1vExMTYf4+Li2t02ic8\nPBxv75+fS5GVlcW4cePsr81mM6GhoWRmZtK3b98LLi8rK4u4uLhG4ywWC2fPnm22ltY0XV79MjMz\nM+2vIyMj7b/7+vpSWlra5uU3nbfhdmlaR8O6G/5+6tQpampqiIqKajjLq0BGg9fZ9b8opcrr9jMz\nEAJ4AGca7HuGJvM2/J26AHkRmAr4101/7kLvszlyZCDarb69IMHFjgrqjesTQEqua99vkJHx8+dJ\neno60dHR9tdNAzA6OppTp07ZX5eVlZGfn0/v3r1bXV7TedPT0zGZTPTq1avZ9bUWvk2XV7/MhrX0\nhKioKE6fPm1/3fD9x8TE4OXlRV5eHoWFhRQWFqKUClBKDW/DojOAKiBMKRVU99N03qZXN/y9bthI\npVQAcDO2U0ftImEg2i05p4Rz5TUkxDp3f0QtGR8T4PL3G7z00kucPn2agoICHn/8cebOndvitDfe\neCMrV67kwIEDVFVVsXjxYi666CL7UQHA008/zblz58jIyODFF1+0L+/GG2/k+eefJy0tjdLSUhYv\nXszcuXMxmZo/KREeHo7BYCA1NbXZ8TNnziQ5OZk1a9ZgsVj44IMPOHLkCNdcc03HN0YH3HDDDbz4\n4otkZmZSWFjIk08+aR8XFRXFjBkzuP/++ykuLsZqtVLXyNvqqRul1BngW+BZTdMCNE0ztGFef6AU\nKNI0rTfw5468JwkD0W47Trhme0G9YZF+eJk0dqUV6F1Kt7npppuYMWMG8fHx9O/fn0ceeaTFaa+4\n4gr+9re/8atf/YqoqChOnDjB+++/32ia6667joSEBMaMGcOsWbO47bbbAPjNb37DggULuPTSS+nX\nrx/e3t7861//anFdvr6+PPzww0yePJmgoCB27NjRaHxoaChffvklzz77LKGhoTz11FN8+eWXhIWF\ndWJrtN8dd9zBjBkzGDVqFGPHjmXmzJmYTCaMRiMAb7/9NtXV1fYrtoCPsbULtMVCwBM4gu10T2vz\nLgPGAUXAWuC/HXlPWjuvp5aLrwW/f3cv+9ML+PK3Y/Uupdss+iCJ4ioL6+6d5hztBitn2f5/69pW\nJ+3bty9vvPEGV1zR6k2pbaJpGsePH2fAgAFdsjxn9PXXX7No0aLzTmE14PA7kRwZiHZRqr69wLUu\nKW0qIcaf4znlFJQ2f525cG8VFRV89dVXWCwWMjMzWbZsGXPmzNG7rE6RMBDtkpJTSn5Ztct0TteS\ncXXtBjtcuN1AdJxSiiVLlhAcHMzYsWMZOnQojz32mN5ldYpcWiraZUfdeXRXu/O4qRFRZjyNGjvS\n8pk1NlbvcrpUw/sDuoI7dt3h6+vL7t279S6jS8mRgWiXnan5hJs96RPkWvcXNOVlMjAiysyeU0Vu\n+WEn3I+EgWiz+vaCcX38naNRtZMSYgI4llNGYZm0G3SFadOm8cYbb3Ro3vT0dMxmM7W1tV1clagn\nYSDa7GR+OTklVS7fXlAvIcYfq4IdJ5rv6Ex0n759+/Ldd9/ZX8fGxlJaWmq/dFN0PQkD0WY76/rr\ncfUrieqNqHsu8k4X76dICJAwEO2wM62AEF8P4kK8W5/YBXh7GBgeaWa3Cz7f4ELdWHe0u+qlS5dy\n880326dt7kE19U6cOMHll19OaGgoYWFhzJ8/n8LCQgAWLFhAeno61157LWazmaeeeuq8ZWVlZTF7\n9mxCQkIYMGAAr7/+un3ZS5cu5YYbbmDhwoX4+/szfPhw9uzZ0/Ub0cVIGIg225mWz1g3aS+olxDj\nz9GzZRSWVeldSpdrrhvrznRX3R5KKR566CGysrJISkoiIyODpUuXAvDOO+8QGxvLF198QWlpKQ88\n8MB588+bN48+ffqQlZXFxx9/zOLFi9mwYYN9/Oeff868efMoLCxk9uzZ53WbLc4nYSDaJKOgnKzC\nSpd7fkFrEmICqFWwK9X12g2a68a6M91Vt8eAAQO48sor8fLyIjw8nPvuu4/Nmze3ad6MjAx+/PFH\nnnzySby9vRkzZgy33347b7/9tn2aKVOmMHPmTIxGIwsWLOCnn35qd43uRsJAtMlO+/0F7tF4XG9U\ntBmTQWO7C4ZBc91ON+0iumF31Rear73Onj3LvHnz6N27NwEBAdx8880tPpGsqaysLEJCQvD3/3lf\nbK0b68rKSnmsZiskDESb7EzNJ9DbRHxY6w8gcSXeHkZGRLlmu0Fz3U53prtqPz8/ysvL7eOys+1d\n9p9n8eLFaJpGYmIixcXFvPvuu42274VORUZHR1NQUEBJSUmjOnq6G2tXI2Eg2mRnWgFj+/hjcKP2\ngnoJMf4czXa9doPmurHuTHfVY8aMYcuWLaSnp1NUVMQTTzzR4rpLSkowm80EBgaSmZnJ008/3Wh8\nr169WuzGOiYmhksuuYSHHnqIyspKDh48yJtvvtmo8Vq0n4SBaNWZogrSC8rd5v6CpurbDXamulY/\nRc11Y92Z7qqvvPJK5s6dy6hRo0hISLjgMwaWLFnCvn37CAwMZNasWVx//fWNxj/00EMsX76coKCg\nZp9h/N5773Hy5Emio6OZM2cOy5Yt67JeWN2VdGEtWvXZgUzufv8A7y4YwZBefnqX0+Mqa2qZvmIv\n88dHs3TOGL3LaV47urCGjndjLd1Vd5jDH1LLkYFo1Y7UAsxeRgaG++pdii5cud1AiHoSBqJVO9Py\nGdPbH6PB4b/cdJuEPgEkuej9BkKAhIFoRU5xJam5ZW53f0FTCbH1/RS5RrvByZMnO3SOXSklp4hc\nlISBuKDtdf3yjHfzMBhZ10/RDumnSLgoCQNxQTtS8zF7GRkU4Z7tBfW8PQyMlHYD4cIkDMQFbT9h\n64/IndsL6iXEBNT1UyTPNxCuR8JAtCi7qJKT+eVuf4qonjzfQLgyCQPRovr+eNy98bjeyGh/PI0a\n209Iu4FwPRIGokU7ThQQ4G1y+/aCevXPRZZ2A+GKJAxEi7an5rttf0QtqX8u8rlSaTcQrkXCQDQr\ns9DWH5G0FzSWEBMg7QbCJUkYiGbtqDsvPj5WwqChkdFmvEwa2yQMhIuRMBDN2p6aT5CPif5u9vyC\n1niZDIzu7c+OtHPSbiBcioSBaFb9/QXSXnC+ibGBpORVcLaovPWJhXASEgbiPBkF5WQWVjAhNlDv\nUhzShLpTZz8k5+hciRBdR8JAnKe+PyJ3e95xWw3p5Ye/l5EfU6TdQLgOCQNxnh0n8gn2lfaClhgN\nGuNjA9h5shCr1ap3OUJ0CQkD0YhSiu2p+ST0CbjgQ8nd3cTYQM4UV5OaW9L6xEI4AQkD0cip/HLO\nFFWSIJeUXtCEuLp2g2PSbiBcg4SBaOSHuvPgEyUMLigu2JsIswfbpJ8i4SIkDEQjW4/nEhngSWyw\nt96lODRN05gYF8iu9CIsllq9yxGi0yQMhJ2l1sq2E/lcFBco7QVtMCE2gKIKC4mnC/QuRYhOkzAQ\ndgcziyiptHBRnNxf0BYT67bTD8mu8Vxk4d4kDITd1uO29oIJ0l7QJuFmT/qFeLM9TY4MhPOTMBB2\nW4/nMTjCl2BfD71LcRoT4gLZn1FCRXWN3qUI0SkSBgKA0ioL+9LPySmidpoYG0ClxcruVLkbWTg3\nCQMBwM7UfCxWxUV9JQzaIyEmAINmuwpLCGcmYSAA+OF4Hl4mA2N6S39E7eHvbWJ4pJkfUwukS2vh\n1CQMBABbU/IY09uMl0l2ifa6pF8gR86UkVNcoXcpQnSY/MsXnCmqICWnlEl9g/QuxSld0i8IBWxM\nyta7FCE6TMJA2C8pnRgnl5R2xNBIP4J9TGyW+w2EE5MwEGxNySPY18TAcF+9S3FKBk3j4n5BbEsr\npEa6phBOSsLAzVmtiq0peUyMDZRHXHbCJf0CKaqwsP+UdFwnnJOEgZs7ml1Cfmm13F/QSZP6BmLQ\nYOPRs3qXIkSHSBi4uS1118df1FfaCzojyMeD4ZFmthzPl0tMhVOSMHBzG47mMCjCl17+XnqX4vQm\nxwdxJFsuMRXOScLAjRWWV7P31DmmxMslpV3hkn6BcompcFoSBm5sc3IutVbFVJ3CIDk5mTlz5hAf\nH8+ECRNYu3YtAOnp6YSHhxMXF2f/efbZZ+3zrVixgsGDBzNlyhSOHDliH75z504WLlzY4++j3pBe\nfoT4yiWmwjmZ9C5A6Of7pByCfU0MizT3+LotFgsLFizglltu4eOPP2bbtm3cfPPNbNiwAQ8PW6+p\nJ06cwGRqvItmZ2ezevVq9uzZwwcffMDy5ctZs2YNFouFJUuW8Nprr/X4e6ln0DQm9Q3ix7pLTD1M\nRt1qEaK95MjATVlqrWw6lsPkfkEYDT1/Senx48fJzs5m0aJFGI1Gpk6dysSJE/nwww8vOF9mZiYj\nR47E39+fyy67jFOnTgHw6quvctVVVxEbG9sT5bdocnyQXGIqnJKEgZvae+ocxZUWpvYP1rsUO6UU\nR48etb8eO3Yso0aN4o9//CP5+bYP1379+pGUlERRURGbN29m8ODBZGZm8sknn/CHP/xBr9LtJsXZ\nLjHdkCSXmArnImHgpjYczcFk0LhIpy4oBgwYQHh4OCtWrKCmpoaNGzeybds2ysvLCQkJYf369ezf\nv5/vvvuO0tJSFi1aBEBISAj33nsvc+bMYf369SxbtoyHH36YRx99lLVr1zJ79mwWLFhAVlaWLu8r\n0MfEiCgzW1LkElPhXCQM3NT3R3MYF+OP2UufZiMPDw/+/e9/s379eoYPH87LL7/MddddR3R0NGaz\nmTFjxmBtnrNoAAAVbElEQVQymYiIiOAf//gHmzZtorS0FIDrr7+eDRs28MEHH5CUlISnpycjR45k\n6dKlrF69mtmzZ7NkyRJd3hfA1P62S0xPF5TpVoMQ7SVh4IZO5ZeRklPK1Hh9TxENHz6czz//nOTk\nZD766CNOnTrFuHHjzptOq+smw2q1NhpeUVHB3//+dx577DFSU1Pp3bs3/v7+jB07ttFVRj1t+sAQ\nAL5OzNStBiHaS8LADW04mgPYvsHq6fDhw1RWVlJeXs5LL73E2bNnmTdvHnv37iUlJQWr1UpBQQGL\nFy9m8uTJBAQ0PqX13HPPMW/ePCIjI+nduzcpKSnk5OSwdetW4uLidHpX0DfEh/hQH749kqNbDUK0\nl1xa6oY2HM2hb4gPfYK8da3jo48+4t1336WmpoZJkybx0Ucf4eXlxalTp3j88cfJy8vDbDYzbdo0\nXn311UbzHj9+nE2bNvHNN98AEBkZyV133cXUqVMJCwvj9ddf1+Mt2U0fGMzKnVnkFlcQHuCjay1C\ntIXWzkYuaRFzcqVVFsY+9i3zxkZy9zR9L8N0ZUfPlnHzO4f427VDWDC5f/evcOUs2/9vXdv96xId\n4fBdAstpIjez9XguNbWKKTqfInJ1gyN8iQ7w4tsjcompcA4SBm5m/ZEc/L2MjI7u+buO3YmmaUwf\nGMyOk4UUllXpXY4QrZIwcCPVFivfHsnm0v7BmIzyp+9u0weFUFOr+O6wPvc8CNEe8ongRn44nktJ\npYUrh4TqXYpbGBVtJtTPg3WH5VSRcHwSBm5k7cEz+HsZdbvr2N0YNI1pA4LZmnqO8qoavcsR4oIk\nDNxEZU0t3x45y7SBIXjIKaIeM31gCBU1VjYmndG7FCEuSD4V3MSW5FxKqyzMGByidyluZXyMP/5e\nRr45JKeKhGOTMHATaxPPEORjYkKsnCLqSSajgUv7B7M5JZ+qGove5QjRIgkDN1BZU8t3R84yfaBc\nRaSHXwwOobiylg1H5HGYwnHJJ4Mb2HQsh7LqWq4YLFcR6eGSvoEE+Zj47/7TepciRIskDNzAlwfP\nEOzrQUKMnCLSg8lo4JdDQ9l0vIBzcgOacFASBi6uvNrC90k5XD4wGJMOj7cUNtcMD6emVvHpvnS9\nSxGiWRIGLm7j0Vwqamq5Uq4i0tXgCF/iQ334ZH+WPAFNOCQJAxf35cEsQv08GNtHThHpSdM0Zg0P\n42BWKSfOFutdjhDnkTBwYefKqvn+aA5XDArBKKeIdDdzWBgGDT7cLaeKhOORMHBh/92fSbXFynUj\nw/UuRQDhZk8mxgXyRWI2tbXW1mcQogdJGLgopRTv70pneJQfgyL89C5H1LlmeBhniqvZliKPxBSO\nRcLARe1LP8fxnFLmjIzQuxTRwLQBwfh5Gvh4b4bepQjRiISBi3pvVwa+nkZmSHfVDsXbw8gvBoWy\nPimP0spqvcsRwk7CwAUVV9bw5cEsrhoSiq+nUe9yRBOzhodRXmPlc7kjWTgQCQMX9Nn+TCprrMwZ\nJQ3HjmhcH3/iQ314e3s6Vqs0JAvHIGHgYpRSvLcrg8ERvgztJQ3HjkjTNG5MiORoThnbU3L1LkcI\nQMLA5SRmFnHkTDH/MzICTZN7CxzV1UPDCPQx8ebWNL1LEQKQMHA57+3KwNtk4Oph0nDsyLw9DFw/\nKoKNyfmk5ZboXY4QEgaupLTKwucHMrlycAhmL5Pe5YhW/HpMBAaDxls/nNC7FCEkDFzJ6h2nKKuu\n5f+N7aV3KaINevl7ccWgEP574AxF5dK1tdCXhIGLqKyp5Y2taVwSH8KwSLPe5Yg2ujEhkrJqK+/v\nPKV3KcLNSRi4iI/3nia3pIrfTo3TuxTRDiOizIyKNvPOznQs0l+R0JGEgQuw1Fp5ZfMJxsYGcVHf\nYL3LEe1047hIThdWse5Qpt6lCDcmYeACvjiYxelzFfxh2gC5nNQJTR8UQqS/Jy9vSpUH3wjdSBg4\nOatV8fLGEwyJ9OfyIdIpnTMyGTR+e0kfDp8p5ZtD2XqXI9yUhIGT+y7pLMdzSvn9tP4Y5AE2Tmvm\n8DD6hfrw9LfHpO1A6ELCwIkppXhp0wliQ3yZNTJK73JEJ5gMGndeGktqbhn/2Scd2ImeJ2HgxLYc\nz+OnjEIWXdYfk1H+lM5u2oBgxsQE8cJ3x6msqdW7HOFm5BPESdXUWvnbl0eIDfHlVwm99S5HdAFN\n03jwl0M4U1TJO9vlvgPRsyQMnNTb20+RklPKo9cMw8skzyxwFRf3D+XSQeG8tCmF4soavcsRbkTC\nwAnlllTxwvpkLhsUzi+GyhVEruaBqwZTWF7Da5tT9S5FuBEJAyf01DdHqbTU8ui1w+S+Ahc0oncg\n142J5rUfUknJkR5NRc+QMHAyBzIK+WjvaX4zuR/9w6UPIlf1yKxh+Hoa+dNHB6m1yo1oovtJGDgR\nq1Wx5LNDhPt78cdfDNS7HNGNwv29WDZ7OAcyCnlzq5wuEt1PwsCJfLgng59OF/HQ1UPkeQVuYPbo\naGYM68Uz3yaTklOqdznCxUkYOImUnFIe+/IIF/UL4X/GyKWk7kDTNJbPGYGPh5EHPv5JTheJbiVh\n4AQqa2q5c80+vD2MvDhvrHQ74UYi/L1ZNns4+9ILWfmjPC9ZdB8JAyew7IvDHM0u4bkbRhMZ6K13\nOaKHXTcmmiuG9uLpdcf4KaNQ73KEi5IwcHCfHcjkvV0Z/H5af6YNlnsK3JGmafzjVyMJ9/fi9rf3\nkFlYoXdJwgVJGDiw1NxSFv83kfFxwdx/5SC9yxE6CjN78dYtE6isruW2VbsprbLoXZJwMRIGDqqg\nrJpF7+7Fw2TgnzeOlY7oBIN6+fPS/HEczynlrvf2S4Oy6FLyCeOAispruPmNnZzKL+flm8YRHeSj\nd0nCQVw6KJyls4ez4WgOy9ce0bsc4ULkYnUHU1JZw8K3dpKSU8prCxO4ZECY3iUJB7NgUhxpuWW8\n9WMafp4m7p8xCLm+THSWhIEDKauycMvK3RzOKuaVmxOkwVi06OFZQymrsrBiYwq5JVX8A4UmkSA6\nQcLAQRSWV/Pbd/ZyIKOQl24ayxXDeuldknBgRoPtCqNeAV78c0MKtwaWMLCXP9KZuegoCQMHcPB0\nIb9/d5+ta+q5Y/jlCHmEpWidpmncN2Mw4QHeFH5VQ9KZYqLLqgnx89S7NOGEpAFZR0op3tlxil//\n33YAPlp0MdeOjta5KuFsFkyKY1CEmbJqC1e9sIX1R87qXZJwQnJkoJOiihqWfn6YT/ZnMm1wOM/f\nMIZg+UYnOijEz4sRHkbCKr244+09XD+2N0uuHU6gr4fepQknIWHQw2pqrby3K53n1ydTWFHDfVcO\n4s7pA6S/IdFpfp4mPrttMis2pvDSxhR+PJHH0muHc9XwSNm/RKskDHqIUoqNx3J4fG0SJ3LLuDg+\nlEeuGcrw6EC9SxMuxNNk4L4rBzFjWC/+9NFP/H71PoZE+vPHywdy9QgJBdEyTal23cUotzy2U2mV\nhc8PZPHernQSM4uID/Nj8cyh/GJoRLc8srKyspLSUun73hkZDAZCQkI6NvPKWbb/37rWPshSa+Xz\nn7JYsTGF1NwyBkSYWXRZf64eEYmfPA+jpzl8CksYdANLrZX9GYX8d18mnx/IpKy6liGR/iy4OI4b\nxsfg0Y1dS0gYOK+uDoN6tVbFV4lnWLEhhWNnS/DxMPKLoRHMHh3NZYPD8TLJBak9wOHDQL4edAGl\nFKfPVbA1JY8tyblsTcmjpNKCt4eBa0dFc+NFsYyNCZKH1wtdGA0a146OZtbIKHafLOCLg1l8lZjN\nlwfP4O9lYmK/EPvPiN6B3fplRTguOTJoh1qrIqekktPnKjiVX07SmWIOZxVxJKuY4kpbL5JRgd5c\nOjCcSweFM3VQGAHePXs1hxwZOK/uOjJoTk2tlR9T8lh3OJudaQWk5pYB4ONhZHCkPwMjzAzsZWZg\nhD9xob5EBfrg4ylHEJ3g8N8EXSoMlFIoBValqFUKqxVqlaLWqrDUWrFYFTW1Viy1iupaK1U1Vqos\ntVRZrJRVWSirtlBaVUtZlYWiihrOlVWTX1ZNQVk1uSVVnCmqoKb2503g7WFgSGQAw6IDGBYVwEX9\nQhgQYdb1CEDCwHn1ZBg0lVtSxe6TBexKK+BYdgnHc0rJK61qNE2gjwdRgd5EBHgT5ONBkK8HQT4e\nBPp64udpxNfLhK+HEV8vI14mI14mA14mA54mAx5GAyaDhtGgYTIYMBo1jJqGptmOXAyahkHDlY+e\nHf6NtSsMhj/6TafCoC0zNy1HoezDVIOF1A9X1IVAM/N2hsmgEeznSaifJyF+noSavegT7FP340tM\nsA9xoX4YHezqDAkD56VnGDTnXFk1KbmlZBSUc6aokuyiSrKLK8kprrR9WSqvobiypkv/3QFoGhi0\nn3ta0jRs/S41+KfWcJzttdbodbvW1+FK2+7wY790rA+KZrQrDAYmTFHlxT332L2KinJ8fHzP/2Np\nF3xp/3ah1f1Hw7azaPW/130DabjTGeq+ndi+pdiGd1Rubi7h4eEdnr8zlFLU1tZ2eP78/HxCQ0O7\nsKKe4ax1w8+1a5qG0djBUzF5x23/DxvYdYW1on4/r7UqrEphVWBt8Hv9lzRr3be2+i9sCmX/Ymf/\n9Kkb3uBlYxf4mOpIFtV/tvSUrJTD65RSv+yxFXaAQ58mGj9+PHv27OnJVXYJZ60bnLd2Z60buqj2\nbjgyaI3bb/P2cfgjA7lsQAghhISBEEIIMC5durQ907dr4q6QkJDQ06vsEs5aNzhv7c5aN3RB7QfW\n2P4/dn7ni2kHt97m7bOsJ1fWEQ7dZiCEaCMd2gxEu0ibgRBCCMcnYSCEEKJ7w6CgoIA5c+bg5+dH\nXFwca9asaXa6559/nvj4eAICAoiOjubee+/FYrF175Ceno7ZbG70o2kazz77LACbNm3CYDA0Gv/v\nf/+7x2qvV11dzdChQ+nTp0+j4QcOHCAhIQFfX18SEhI4cOCAfZxSigcffJDQ0FBCQ0N58MEHaedp\nu26pOzk5meuuu47w8HBCQkK46qqrOHbsmH38qlWrMBqNjbb5pk2bOlV3V9UOtntI/Pz87LXdfvvt\n9nGOus1/+OGHZvfz//znP4D+23zp0qV4eHg0Wn9qaqp9vKPu5xeqW6/93GHZunBo80+7zJs3T91w\nww2qpKRE/fDDDyogIEAdOnTovOlSUlLUuXPnlFJK5efnq+nTp6tnn3222WWmpqYqg8Gg0tLSlFJK\nbdy4UfXu3bu9pXVZ7fWWL1+upk6d2qiWqqoqFRsbq5577jlVWVmpXnzxRRUbG6uqqqqUUkq98sor\natCgQSojI0OdPn1aDR06VP3f//2f7nXv3LlTvfHGGyo/P19VV1erRx55RA0ePNg+fuXKlWry5Mmd\nqrO7aldKKUAdP3682XkcdZs3tXHjRmU2m1VpaalSqg3b/K2Ztp9uqn3JkiVq/vz5zS7DkffzC9Xd\nw/t5ez9re/yn+xYMfkA1MKjBsHeAf7QyXyjwHfByC+OXABsbvJ4GnNazdqAfkARc3bAWYAaQSV1D\nfd2wdOCXdb9vA37bYNxtwA69625muhBsFw+E1r2+BdjqiNu8bpwCBrQwn7Ns85XAygavdd3m2K4k\nfLeF5Tjsfn6hupuZttv3c0f+6c7TRIMAi1IqucGwn4DhzU2sadpNmqYVA3nAaODVZqbRgIVA0/NA\nEZqmndU0LU3TtOc1TfPrydqBfwGLgYomw4cDB1XdnlXnYIPlDK9bblvW0RZdVXdTlwLZSqn8BsPG\napqWp2lasqZpf9U0rbPdoXd17Vs0TcvWNO2/mqb1bTDc4bd53f77a87fz/Xe5tdqmlagadphTdN+\n32C4o+/nLdXdVE/s5w6rO8PADBQ3GVYE+Dc3sVJqjVIqANsf+hXgbDOTTQF6AR83GHYUGANEAZcD\nCcBznaq8HbVrmjYHMCqlPmlhOUUXWE7T8UWAWet4141dVXfD6foALwH3NRi8BRgBRAC/Am4E/tzB\nmut1Ze2XAX2BIUAW8GWDf8QOv82B67F9KdrcYJiu2xz4EBgKhAN3AI9qmnZjg+U45H7eSt12Pbif\nO6zuDINSIKDJsACg5EIzKaWOA4eBl5sZ/b/Af5RSpQ2mz1ZKHVFKWZVSacAD2P5wndGm2uu+wT0F\n3NXB5TQdHwCUNvmG1R5dVXf9dOHAt9hO2b1XP1wplaqUSqvb5onAY9i+yXZGl9WulNqilKpWShUC\nd2M7NTO0hfU41Dav87/A2w1r0nOb163/iFIqSylVq5TaBrzYYP0OuZ+3oW6gx/dzh9WdYZAMmDRN\na9iN4mhsH/StMQH9Gw7QNM0H+H+cf+jclKLz76uttQ/E9g30B03TsoH/AlF1pyf61k0/qsk3oFEN\nlnO4brkXWocedaNpWjC2fyCfK6Ueb2W9is7fVNNltbdSn8NucwBN02KwtYO93cp6e3Kbt7Z+R93P\nm9Nou+mwnzuu7myQAN4H3sPW4DMZ26Hc8Gamux2IqPt9GLY/6nNNprkJOEmDRqq64dOBOGx/pBhg\nIw0a3rqzdmyhFdng53pspyUiASPgCZzC9u3UC7iz7rVn3fyLsDUm9gai6973IgeoOwDYBaxoYR1X\nA73qfh8CHAKWOMg2H47ttKER2+mEF4BjgIcjb/MG0y0GtjjSNq+b7joguO7f2URsDcb/WzfOIffz\nNtSty37uqD/du3Bb6/ynQBm2qwtuqhs+FdthYv10K7G1EZRh+8B/GvBusqx1wN+aWcd9dX/gciAD\n+Cfg31O1N5lnGudf2TIW2IutwXAfMLbBOA3bqYOCup+naBJ2etSN7TSFqltGaYOf2LrxzzT4e6Vi\nO3z2cIRtjq3d6FjdMnLqljfQ0bd5g+FHgduaGa7rNsf2wZtftx8cBe5yhv38QnXrtZ876k97+yYS\nQgjhgqQ7CiGEEBIGQgghJAyEEEIgYSCEEAIJAyGEEEgYCCGEQMJACCEEEgbCBWiaNlbTtB81TSvX\nNG2XpmmxetckhLORMBBOra63ya+AJ7E9CyMVeETXooRwQhIGwtk9C7yulPpcKVWBrc+aCTrXJITT\ncdkHNQjXp2laALaOyAY1GGwAKvWpSAjnJWEgnNkvAA/gYIPek72Az3SrSAgnJaeJhDPri60f+qD6\nH2xdmH+jb1lCOB8JA+HMvLB1XQ6Apmn9gPHA57pVJISTkjAQzmw3cJmmadF1TwhbAzyslCrQuS4h\nnI60GQhntgH4EttjEPOBJ5VSr+tbkhDOSR5uI4QQQk4TCSGEkDAQQgiBhIEQQggkDIQQQiBhIIQQ\nAgkDIYQQSBgIIYRAwkAIIQQSBkIIIYD/DxkSYou68UE8AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1b03d89518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create grid of 80 points from 0.36 to 0.54\n",
    "x = np.linspace(0.36, 0.54, 80)\n",
    "# freeze a beta distribution object with given parameters\n",
    "dist = beta(438, 544)\n",
    "# probability density function at x\n",
    "pd = dist.pdf(x)\n",
    "\n",
    "print('Uniform prior -> Posterior is Beta(438,544)')\n",
    "\n",
    "# plot pd\n",
    "plt.plot(x, pd)\n",
    "# show only x-axis\n",
    "plot_tools.modify_axes.only_x()\n",
    "\n",
    "# annotate the line\n",
    "plt.annotate(\n",
    "    r'$p(\\theta|y,n)$',\n",
    "    (x[35] - 0.005, pd[35]),\n",
    "    ha='right'  # horizontalalignment\n",
    ")\n",
    "\n",
    "# plot proportion of girl babies in general population as a vertical line\n",
    "# ``color='C1'`` corresponds to default color #2\n",
    "plt.axvline(0.485, color='C1')\n",
    "# annotate the line\n",
    "plt.annotate(\n",
    "    'proportion in general\\npopulation',\n",
    "    (0.485 + 0.005, 14),\n",
    "    ha='left'  # horizontalalignment\n",
    ")\n",
    "\n",
    "# find the points in x that are between 2.5% and 97.5% quantile\n",
    "# dist.ppf is percent point function (inverse of cdf)\n",
    "x_95_idx = (x > dist.ppf(0.025)) & (x < dist.ppf(0.975))\n",
    "# shade the 95% central posterior interval\n",
    "plt.fill_between(x[x_95_idx], pd[x_95_idx], color='0.92')\n",
    "# add text into the shaded area\n",
    "plt.text(dist.median(), 8, \"95%\", horizontalalignment='center')\n",
    "# add labels and title\n",
    "plt.xlabel(r'$\\theta$')\n",
    "\n",
    "# scale x-axis tightly to the data.\n",
    "plt.autoscale(axis='x', tight=True);\n",
    "# N.B. the last semicolon is here just to prevent ipython notebook\n",
    "# from displaying the return value of the last command."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
